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ABSTRACT 

We use local numerical simulations to investigate the strength and nature of magnetohydrodynamic 
(MHD) turbulence in the outer regions of protoplanetary disks, where ambipolar diffusion is the domi- 
nant non-ideal MHD effect. The simulations include vertical stratification and assume zero net vertical 
magnetic flux. We employ a super time-stepping technique to ameliorate the Courant restriction on 
the diffusive time step. We find that in idealized stratified simulations, with a spatially constant am- 
bipolar Elsasser number Am, turbulence driven by the magnetorotational instability (MRI) behaves 
in a similar manner as in prior unstratified calculations. Turbulence dies away for Am < 1, and 
becomes progressively more vigorous as ambipolar diffusion is decreased. Near-ideal MHD behavior 
is recovered for Am > 10'^. In the intermediate regime (10 < Am < 10^) ambipolar diffusion leads 
to substantial increases in both the period of the MRI dynamo cycle and the characteristic scales 
of magnetic field structures. To quantify the impact of ambipolar physics on disk accretion, we run 
simulations at 30 AU and 100 AU that include a vertical Am profile derived from far ultraviolet (FUV) 
ionized disk models. These models develop a vertically layered structure analogous to the Ohmic dead 
zone that is present at smaller radii. We find that, although the levels of surface turbulence can be 
strong (and consistent with constraints on turbulent line widths at these radii), the inferred accretion 
rates are at least an order of magnitude smaller than those observed in T Tauri stars. We speculate 
that this discrepancy may be due to the neglect of vertical magnetic fields. If this is the case, then 
the MRI alone may result in disjoint classes of disk evolution, with those disks lacking net fields being 
weakly viscous at large radii. 

Subject headings: accretion, accretion disks — (magnetohydrodynamics:) MHD — turbulence — 
protoplanetary disks 



1. INTRODUCTION 

The structure and evolution of protoplanetary disks 
play a crucial role in the formation of stars and their 
planetary systems. Disk gas is observed to accrete onto 
the central star at rates that require some form of an- 
gular momentum transport substantially stronger than 
that provided by molecular viscosity. Turbulence has 
long been suggested as the sour ce of enhanced trans- 
port (jShakura &: SyunyaevI 119731 ). This turbulence not 
only allows accretion, but can also play an important 
role in the formation and subsequent evolution of plan- 
ets. At early times, turbulence could act to inhibit 
dust settling and largely determine the coUisional ve- 
locities that affect the balance between fragmentation 
and coagulation of t hese p articles (Ornicl & Cuzzi 2007; 
lYoudin fc Lithwickl l2007t iBirnstiel et al.l l201Il) . Per- 
sistent pressure maxima predicted by some turbulence 
models (iBarge fc Sommerialll995t iJohansen et al.l [20091 : 
lUribe et al.ll201lHSimon et al.ll2012D" may act to concen- 
trate particles, enhancing their coagulation into larger 
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particles. Once planetesimals have formed, gravitational 
coupling to turbu lent fluctu ations in the disk may af- 
fect their growth (|Ida et al.l i2008'). Finally, the strength 
and nature of turbulence determines whether the critical 
co-orbital contributio n to the Type I migratio n torque 
remains unsaturated (IPaardeko oper et al.ll201l"l ). 

At a minimum, turbulence in protoplanetary disks will 
be generated in regions where the m agnetorational insta- 
bility (MRI; [BaMsAH^ie3|l99^) operates. Indeed, 
simulations of the non-linear evolution of the MRI un- 
der ideal magnetohydrodynamic (MHD) conditions yield 
sustained turbulence that transports angular momentum 
outward a t rates in g eneral agreement with observations 
( O^artmann et al.|[T998.) . However, large regions of proto- 
planetary disk s are expected to have very low ionization 
fractions fe.g.. Illgner fc NelsonI 120061 ) . which in turn re- 
sult in three significant non-ideal MHD effects: Ohmic 
diffu sion, a mbipolar diffusion, and the Hall effect (see, 
e.g., lArmit agc 20II). The relative importance of these 
effects depends primarily upon the density (as well as 
magnetic field strength). Ohmic diffusion is efficient at 
high densities, and is thus most important in the inner 
regions of the disk (outside a small zone very close to 
the star where thermal ionization of alkali metals pro- 
vides sufficient ionization throughout the disk column). 
At low gas densities, such as in the outer disk, ambipolar 
diffusion becomes dominant, while at in termediate den- 
sities , the Hall term is important fe.g.. iKunz fc BalbusI 
[2001 . 

MRI physics and the phenomenological consequences 
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of the non-ideal terms have been best-characterized in 
the case of Ohmic diffusion. The evolution in this 
limit depends upon the Elsasser number, defined as 
A = v\^/riil, where va,z is the Alfven velocity in the 
vertical direction, rj is the Ohmic resistivity, and D, is the 
angular frequency of Keplerian rotation. For A less than 
order unity, MRI turbu lence is severely quenched pini 
ll996HTurner et al.|[2007[) . while for larger values the MRI 
saturation level will depend on the strength of Ohmic dif- 
fusion; stronger diffusion leads to lower turbulence lev- 
els. Combining these results with chemical models for 
disks motivates the dead zone model of disk accretion 
(|Gammid 119961 ). In this model, the disk is well- ionized 
only in its surface layers due to non-thermal sources (X- 
rays, cosmic rays, and far ultraviolet (FUV) photons) 
that penetrate the disk from the exterior down to some 
column depth. Closer to the mid-plane, the ionization 
fraction is low, re sulting in a smah A, and no MRI-driven 
turbulence (e g.. lGamini3ll996t iFleming fc Sto^ 120031: 
Turner et al.) I2007t iTurner fc Sand l2008t lOishi fc Lo^ 
20091) . Thus, MRI-driven accretion occurs in active lay- 
ers only, leaving much of the disk mass near the mid- 
plane magnetically inactive. 

Qualitative changes to the predicted disk structure 
may equally result from ambipolar diffusion and the Hall 
effect, though for these terms the understanding of the 
non-linear behavior is incomplete. The linear regime of 
the MRI in the presence of the Hall term has bee n ex- 
plored byjW ardlc ( 1999l .galbus fc Terque^ (pOOl . and 
IWardle fc Salmeron (201^ . A primary result from these 
studies is that the growth rate of the MRI is strongly 
affected by the sign of $7 • B, i.e., how the vertical mag- 
netic field is aligned with the angular velocity vector. 
The only study of the non-linear, turbulent state of the 
MRI in the presence of th e Hall term wa s carried out in 
ISano fc Ston"a (|2002ah and HanoXStoni (|2002H) . Their 
numerical simulations included both Ohmic diffusion and 
the Hall term, with the Ohmic contribution dominating 
significantly over the Hall effect. In this regime, the Hall 
term does not strongly influence the saturated state of 
the MRI. However, the regime in which the Hall term 
dominates has yet to be explored through simulations 
(IWardle fc SalmeronI 12011) . 

Ambipolar diffusion arises from the imperfect coupling 
betwee n ionized species and ne utrals. The line a r ana l- 
yses of Blacs fc Balbus (11994 . iKunz fc Balbui (|200l . 
and iDesch (.20043 showed that the growth of the MRI is 
damped when the collision frequency between the neu- 
trals and the ions is smaller than the orbital frequency. 
This is intuitive; neutrals need to communicate with the 
ions faster than the timescale over which the MRI acts 
(i.e., the dynamical one) in order for the neutrals to feel 
any MRI-like effect at all. 

Initial two and three-dimensional simulations of non- 
linear MRI turbulence in the presence o f amb ipolar 
diffusion were car r ied o ut by iLow et all ()1995l ) and 
[Brandenburg et aTl ()1995[) . respectively. These authors 
considered the single-fiuid, "strong-coupling" limit, valid 
when the recombination timescale is much shorter than 
the dynamical time. T his limit is general l y applicabl e 
to protoplanetary disks (|Bai fc Stonell2011l :' lBail(2011af ). 
Thei r results agreed with t he expectations of linear the- 
ory (jBlaes fc BalbusI |1994[ ): the MRI only operates if 



the collision frequency between neutrals and ions ex- 
ceeds the angular frequency. Simulations in the alter- 
nate regime, where the recombination timescale for elec- 
trons is assumed to b e very long, were conducted by 
IHawlev fc Stio^ using a two-fiuid approach in 

which the ions and neutrals were evolved separately, only 
interacting through collisions. The primary result of this 
work was that for a collision frequency less than O.Olfi, 
the ions and neutrals behave independently, but for a 
frequency larger than lOOO, the gas behaves as though 
its fully ionized. When both the collision and orbital 
frequencies are comparable, the saturation level of the 
turbulence is primarily controlled by the ion density. 

The importance of ambipolar diffusion for outer pro- 
toplanetary disks, the advent of well-resolved mm-wave 
observations of this region, and advances in numerical 
techniques, all motivate more detailed studies of how am- 
bipolar diffu s ion aff ects the saturated state of the MRI. 
iBai fc Stond (|2011[ ) carried out shearing box simulations 
in the strong-coupling limit to determine how the MRI 
saturation level correlates with dimensionless number. 
Am, defined as the frequency for the neutrals to collide 
with the ions divided by the orbital frequency (see § 12. ip . 
They found that for simulations with a net toroidal mag- 
netic flux, but no vertical magnetic fiux, the turbulence 
dies away for Am < 1, in line with previous studies. 
For sustained turbulence runs, the saturated turbulent 
stresses increase with increasing Am, eventually asymp- 
toting towards the ideal MHD level. However, in the 
presence of a net vertical magnetic fiux, turbulence can 
always be sustained even for Am < 1, assuming that the 
background vertical magnetic flux is weak enough. For 
low Am values, however, the resulting turbulence levels 
are fairly small. 

Following the above lines of investigation, we study 
the effect of ambipolar diffusion on the MRI in the 
outer region of protoplanetary disks by performing nu- 
merical simulations using a more realistic disk structure 
than attempted previously. We include vertical strat- 
ification (absent in all p rior work except for that of 
iBrandenburg et aD ()1995[ )). which has been shown in 
some previous MRI calculations to lead to significant 
qualita tive changes in th e non-linear evolution. For ex- 
ample, iDavis et al.l ()2010[ ) showed that in the ideal MHD 
case, the turbulence properties converge with numer- 
ical resolution, while in unstratified simulations (and 
for a vertical domain size of one scale height or less, 
(Stone, private communication)), the stress level de- 
creases with resolution with no signs of convergence 
( Fromang fc Papaloizoul [20071 ). Similarlv. iSimon et al.l 
(|2011bD showed that in the presence of Ohmic resistivity, 
vertical gravity can lead to large amplitude fiuctuations 
in the stress levels, a behavior that is absent without 
vertical gravity. 

We also aim to translate our idealized understanding of 
the ambipolar-dominated MRI into predictions for tur- 
bulence and accretion in the outer regions of protoplane- 
tary disks. We run simulations that mimic realistic con- 
ditions in the outer disk (which reflect the chemistry cal- 
cu lations in[Ba i (2011a) an d the F UV ionization model 
of iPerez-Becker fc Chiang ()2011bO ). These simulations 
will provide a quantitative measure of the turbulence, 
structure, and evolution of the outer disk regions. 

Our investigation is divided into two sets of studies. 
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using different magnetic configurations. Tlie first, which 
we pursue here (Paper I), assumes zero net vertical mag- 
netic flrufl. Although it is not the most geometrically 
favor able configuration for the ambipolar-dominated 
MRI (|Bai fc Stonell2011[ ). this is the most studied field 
configuration in the li terature (e.g., iStone et"al1 119961 : 
[Fleming fc Stondl2003f ). and it allows us to make direct 
comparison s between the amb ipolar MRI and the ideal 
MHD case (iDavis et al.ll2010D. as we ll as the case with 
Ohmic resistivity feimon et al.ll2011bf ) . It is also entirely 
conceivable (though, maybe not particularly likely) that 
there will be some regions of protoplanetary disks that 
have negligible vertical magnetic fields; our results will 
apply to such regions. The second set of studies will in- 
clude a non-zero vertical net magnetic flux, and we defer 
that problem to Paper II. 

The structure of the paper is as follows. In Section [2l 
we describe our equations, the methods used to solve 
them and the initial conditions for our simulations. In 
Section 13.11 we study the effect of Am (here assumed 
constant in space and time) on vertically stratified MRI 
turbulence. Then, in Section 13. 2[ we apply a realistic 
protoplanetary disk model to allow for a spatially and 
temporally varying Am. Section |4] discusses the implica- 
tions of our results for real protoplanetary disks, and we 
wrap up with conclusions in Section [5l 



2. METHOD 

2.1. Numerical Method 

In this study, we use Athena, a second-order accurate 
Godunov fiux-conservative code for solving the equations 
of MHD. Athena uses the dimensio nally unsplit c orner 
transport upwind (CTU) method of lCoienal ()1990D cou- 
pled with the thir d-order in space piece wise p arabolic 
method (PPM) of iColella fc Woodwar d! (ITOSl an d a 
constrained transport (CT; lEvans fc HawlevI Il988f ) al- 
gorithm for preserving the V • B = constraint. We 
use the HLL D Riemann solver to calculate the nu - 
merical fluxes (|Mivoshi fc Kusano|[2005l: lMignonell2007D . 
A detailed description of the base Athena algorithm 
and the results o f various test problems are giv en in 
Gardiner fc Stond ()2005l ). iGardiner fc Stond (|2008f ). and 
Stone et all ([200g). 

Our set up is specialized and necessarily more complex 
than the base algorithm. First, our simulations utilize 
the shearing box approximation, which is a model for 
a local, co-rotating disk patch whose size is small com- 
pared to the radial distance from the central object, i?o. 
This allows the construction of a local Cartesian frame 
(x, y, z) that is deflned in terms of the disk's cylindrical 
co-ordinates {R,(j),z') via x = {R — Rq), y = Ro<t', and 



z — z'. The local patch co-rotates with an angular veloc- 
ity fl corresponding to the orbital fr e quenc y at Rq, the 
center of the box; see Hawley et al.l ()1995[ ). Thus, the 
equations to solve are: 



dp 
dt 



V • (pv) = 0, 



(1) 



^ Of course due to the dynamo action of the MRI IIDavis et al.l 
120101 : [5imon et al.ll20irED . the net radial and toroidal fields are 
allowed to evolve in time. 



^ + W ■ {pvv - BB) + W (^P +Ib'^ 
= 2qpQ^x — pfl^z — 2Q X pv 



dB 

'dt 



-Vx(t>xJB) = Vx 



{J xB)x B 



(2) 



(3) 



where p is the mass density, pv is the momentum den- 
sity, B is the magnetic held, P is the gas pressure, and 
q is the shear parameter, deflned as q — — dlnil/dlni?. 
We use q = 3/2, appropriate for a Keplerian disk. We 
assume an isothermal equation of state P — pc^, where 
Cs is the isothermal sound speed. From left to right, the 
source terms in equation ([5]) correspond to radial tidal 
forces (gravity and centrifugal), vertical gravity, and the 
Coriolis force. The source term in equation Q is the 
effect of ambipolar diffusion on the magnetic field evolu- 
tion, where pi is the ion density, and 7 is the coefficient of 
momentum transfer for ion-neutral collisions. Note that 
our system of units has the magnetic permeability /i = 1, 
and the current density is 



J = V X B. 



(4) 



Adopting this shearing box approximation allows for 
better resolution of small scales within the disk, at 
the expense of excluding global effects (those of scale 
^ Ro). Thes e scales could be physi cally important 
(jSorathia et al.l 120121: [Simon et al.ll2012D . However, the 
trade-off is worthwhile for our purposes, because we need 
to study not only models where ambipolar diffusion is 
dominant, but also situations where diffusion is only im- 
portant on small scales. 

The numerical integration of the shearing box equa- 
tions require additions to the Athena algorithm, t he de - 
tails of which can be fo und in IStone fc Gardiner! (poTol ) 
and the Appendix of iSimon et al.l i2011b). Briefly, 
Crank- Nicholson differencing is used to conserve epicyclic 
motion exactly and orb ital advection to subtrac t off the 
background shear flow ([Stone fc Gardineill2010D . The y 
boundary conditions are strictly periodic, whereas the 
X boundaries are shearing periodic (Hawl ev et al. 1995). 
The vertical bound aries are the o utflow boundary con- 
ditions described in ISimon et al.l (|2011bl ) . The electro- 
motive forces (EMFs) at radial boundaries are properly 
remapped which guarantees that the net vertical mag- 
netic &miAs_stin£ti^^ to machine precision using 
CT (IStone fc GardinedlMoh . 

The integration of the ambipolar diffusion term also re- 
quires some modiflcations to the algorithm. Ambipolar 
diffusion is ir nplemented i n a flr st-order operator-split 
manner as m iBai fc Ston^ ()201lD : the ambipolar diffu- 
sion term is integrated separately from the ideal MHD 
integrator. Furthermore, as is evident from equation (|3]), 
the ambipolar diffusion term can be written as an EMF. 
Thus integrating it is done via the CT method to pre- 
serve V ■ B — 0. The ambipolar diffusion EMFs are also 
remapped at the radial boundaries in the same way as the 
ideal MHD EMFs in order to conserve the vertical mag- 
netic net flux. In addition, before even remapping these 
ambipolar diffusion EMFs at the radial boundaries, we 
also remap the toroidal current densities Jy (located at 
cell edges) so that the line integral J Jydy along the in- 
ner and outer radial boundaries are equal. We flnd that 
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this procedure is essential to avoid spurious numerical 
features at shearing-box boundariesQ 

Throughout this paper, the strength of ambipolar dif- 
fusion will characterized by the ambipolar Elsasser num- 
ber 

A,n.f. (5, 

which corresponds to the number of times a neutral 
molecule collides with the ions in a dynamical time 
(il^^). Am can be rewritten as 



Am = 



(6) 



which is a form reminiscent of the Ohmic Elsasser num- 
ber. As shown by equation ([SJ, Am is independent of 
the Alfven speed; this comes about because the ambipo- 
lar diffusivity, ?7a is defined as 



VA 



(7) 



This diffusivity is responsible for determining the diffu- 
sive time step in a Courant limited calculation; Atdiflf oc 
Ax^/yyA. Since this diffusivity is proportional to the 
Alfven speed squared, it can become very large in the up- 
per disk regions, making the Courant limited time step 
extremely small in some of our calculations. 

To circumvent this issue, we have i mplemented the 
super time-stepping (STS) technique of lAlexiades et al.1 
(|1996D to accelerate our calculations. The STS technique 
has already been successfully implemented and tested 
for studying ambipolar diffusion in multi-fluid co des by 
lO'Sullivan fc Downesi 12006) and O' Sullivan & Downp 
()2Q07| ) and in a single-fluid code by Choi et al, ()2009( ). 
Our implementation is similar to theirs, as we briefly de- 
scribe below. 

The STS method divides a compound timcstcp Aigxs 
into N unequal substeps At, (j = 1, A'^) with 



N 



At 



STS 



(8) 



By choosing At, judiciously, stability can be maintained 
even when the averaged timestep A^sts /N is much larger 
than the normal stable diffusion timestep A^diff- The 
optimized lengths for the substeps were found to be 
(|Alexiades et al.lll996D 



Ar,- = Atdiff' 



— 1) cos 



2j 



N 



1 TT 

~2 



(9) 



where < < 1 is a free parameter. The sum of the 
substeps gives 



AisTS — Atdiff 



N 



2N 



(1 + + {i-V^) 



2N 



(10) 



* We note that strict m agnetic flux c o nserva tion (remap of the 
EMFs) was not enforced in lBai &: Stone! 120111) . in which case this 
additional remap of Jy was not neces sary. Nevertheless, the vari- 
ations in vertical net magnetic flux in lBai fc Stone! 120111 ') simula- 
tions were tiny (less than 0.01%), which did not affect their results. 



We note that as — 0, AtsTS — > A^^A^diff so that 
the STS approach is asymptotically N times faster than 
the standard explicit approach. However, the value of v 
needs to be properly chosen to achieve the optimal bal- 
ance between performance and accuracy. In general, the 
STS method provides better accuracy as N decreases and 
V increases, whereas large N and small v lead to higher 
efficiency. Here, we choose u = 1/47V^ with a limit of 
A^<12. AtA^ = 12, one achieves an acceleration fac- 
tor of about 9. It is also found that further increasing 
N would not significantly increase the efficiency without 
sacrificing accuracy (Stone, private communication based 
on a Princeton Junior Project done by Sara Wellons). 

In our calculations, we first compute the ideal MHD 
time step AtMHD and the diffusion timestep Atdiff- The 
number of super time steps N can be found from the 
condition G[N ~ 1,1/4(A^- 1)^] < AtMHD/Atdiff < 
G{N,l/m^). li N < 12, then we modify AidifT so 
that AtMHD = G{N,l/AN'^)AUis. Otherwise, we fix 
N ^ 12, and set AiMHD = Atdiff G[12, 1/(4 x 12^)]. 
In this way, we always have AtMHD = AtsTS- As 
we use an operator-split algorithm for magnetic diffu- 
sion, we integrate N STS substeps of the ambipolar dif- 
fusion term with Ar^ before evolving one MHD time 
step with AtsTS- With STS, we have repeated the 
test problems (i.e., linear wa ve damping test an d C- 
type shock test) performed in ' Bai k, Stond (|201ir ) and 
found essentially the same results for A^sts up to 10. 
Moreover, we repeated the unstrati fied MRI simulation s 
with Am = 1 (runs Z5 and M5 in iBai fc Sto^ (|2011[ )') 
and with STS turned on. In these simulations A^sts 
reaches 12, and the stres s level we find is the same as 
reported in IBai fc Stond ( 201 If). Combining our tests 
with the successful tests oflO'Sullivan fc Down cs (2006), 
lO'Sullivan fc Downesi (120071 ) and lChoi et al.l ()2009l we 
are confident that the STS technique implemented here 
is capable of achieving substantial speed-up while main- 
taining accuracy. 

2.2. Am Profiles 

For most of our simulations, we fix Am to be constant 
in order to study the effect of ambipolar diffusion on the 
non- linear saturation of the MRI in the presence of verti- 
cal stratification. This is desi gned to be the next logical 
step in extending the work of IBai fc Stond ()2011f ) where 
vertical stratification was not included. We have con- 
sidered Am=l, 10, 100, 300, 10^ and 10''. Affhough this 
prescription of a constant Am profile is highly simpli- 
fied, it is a necessary, incremental step bet ween the con- 
stant Am models without vertical gravity ()Bai fc Stond 
1201 If ) and simulations that incorporate a more realistic 
prescription for ambipolar diffusion, which we also carry 
out (see below) 

The results of lBail (|2011a[ ) show that the physical value 
of Am in the outer regions of PPDs is o f order unity. Re- 
cently, iPerez-Becker fc ChiangI (j2011bf ) pointed out that 
the surface layer of protoplanetary disks should be much 
better ionized due to far ultraviolet (FUV) photon ion- 
ization from the central star; these photons almost fully 
ionize the carbon and sulfur to overcome the effects of 
recombination onto dust grains. Their results suggest a 
large ionization fraction (/ ^ 10~^) down to a small pen- 
etration depth of Spuv ~ 0.01 — 0.1 g cm~^, relatively 
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independent of disk radius. Such an ionization fraction 
should significantly reduce the strength of ambipolar dif- 
fusion (i.e., increase Am) in the disk surface layers. 

Thus, in our second set of simulations, we include the 
effect of FUV ionization at the disk surface layers to give 
Am a more realistic spatial and temporal de pendence. 
Since we are not including Ohmic resistivity (jGammid 
[19961) in our calculations, these particular models are 
only appropriate for the outer regions of protoplanetary 
disks (e.g., beyond ~ 10 AU). We adopt a minimum- 
mass solar nebular disk model with surface density of 
1700i?Xu^g cm-2 (IWeidenschiUind [19771 iHavashil 



I1981D . where i?AU is the disk radius measured in AU 
We can express the value of Am wi thin the FUV ionized 
layer as follows (|Bai fc Stonell2012[ ) 



Ampuv « 3.3 X 10^ 



/ 



10- 



PO.mid 



R 



-5/4 
AU 



(11) 



where / is the ionization fraction and po,mid is the mid- 
plane density. For simplicity, we fix / = 10~^, po.mid = 1, 
and assume a penetration depth of Spuv = O.lg cm~^ 
(whic h is slightly different from that in iBai fc Stond 
|2012| ). We conduct two simulations that correspond to 
radial locations at i? = 30 AU and R = 100 AU. Assum- 
ing the density profile is Gaussian (see Equation [15]), one 
finds that the base of the FUV layer (at which the col- 
umn density equals Spuv) is located at z^, — 1.1 H for 
i? = 30 AU and = I.IH for R = 100 AU {H is the 
vertical scale height as defined in equation (|16|) below). 
In our simulations, we set Am = 1 for — z;, < z < Zf, as a 
proxy based on the calculations of iBail (|2011aD . and use 
Equation (ITU) for the ionized surface layers the disk. For 
simplicity, we keep the value of zj, fixed throughout the 
calculation. 

From these considerations, the value of Am changes 
quite dramatically from Am=l to Am« 3 x 10^ at the 
base of the FUV layer. This very large transition is 
smoothed over roughly 7 grid zones so as to prevent a 
discontinuous transition in Am. The smoothing func- 
tions we apply are based upon the error function (ERF). 
Thus, the complete profile of Am for these runs is given 

by 



Am = < 



z > Zb + /S.Z 

Zb — n^z < z < Zb + Az 
—Zb + nAz < z < Zb — nAz 
—Zb — Az < z < —Zb + nAz 
z < —Zb — Az 

(12) 

where S~^{z) and S (z) are the smoothing functions de- 
fined as 



Ampuv 

1 + iAmFuv5'+(2:) 
1 

1 + |AmFuvS'"(2) 
Ampuv 



S+{z) = 1-hERF 



S~{z) = 1 - ERF 



z - 0.9zb 
A'z 

z + 0.9zb 
a'z 



(13) 



(14) 



Here, n = 8 and Az = 0.05H. These numbers were 
chosen to give a reasonably resolved transition region 
between Am = 1 and Ampuv- For a visual representation 
of the rather complex equation ([T2|) . we plot in Fig. [T] 
the vertical profile of Am (averaged over x and y) for 



E 
< 




Figure 1. Vertical profile for Am at -Rau = 30. The profile 
corresponds to the initial time of Z30AU, which is orbit 22 from 
the Am = 10^ run with the same domain size. The value of Am 
has been averaged horizontally. The units of the x axis are the 
vertical scale height, H. The asterisks denote the locations of grid 
zones. Am transitions from Am = 1 to Am = Ampuv (^s defined 
in the text), and this transition is smoothed over roughly 7 grid 
zones using the error function. 

the run at i?AU = 30 at the initial time, referring here 
to when the run was restarted from Cle5 (see below). 
The asterisks denote the grid cell centers; as previously 
mentioned, the transition region is resolved by '--^ 7 zones. 

2.3. Simulations 

We have run simulations with several domain sizes and 
Am profiles. All of the simulations with Am < 10^ or 
with spatially and temporally varying Am are initialized 
from the turbulent state of a "starter" calculation with 
the same domain size but with Am = 10^ (i.e., reasonably 
close to ideal MHD). 

These starter simulations are initialized with a density 
corresponding to isothermal hydrostatic equilibrium. 



p{x,y,z) = poexp - 



(15) 



where po = 1 is the mid-plane density, and H is the scale 
height in the disk. 



H = 



V2cs 

n 



(16) 



The isothermal sound speed, Cg = 7.07 x 10 corre- 
sponding to an initial value for the mid-plane gas pres- 
sure of Po = 5 X lO^'^. With n = 0.001, the value for 
the scale height is H — 1. A density fioor of 10~* is 
applied to the physical domain as too small a density 
leads to a large Alfven speed and a very small time step. 
Furthermore, numerical errors make it difficult to evolve 
regions of very small plasma f3 (ratio of thermal pressure 
to magnetic pressure). 

The initial magnetic field is purely toroidal and has a 
constant f3 — 100 throughout the domain (thus, By has a 
Gaussian shape like the density) . Random perturbations 
are added to the density and velocity components to seed 
the MRI. 
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The remaining simulations are restarted from orbit 50 
(orbit 22 for the simulations with domain size 8Hx 16Hx 
8ii43) of their corresponding domain size simulation with 
Am = 10^ They are hsted in Tabled The label for 
each calculation describes whether the value of Am is 
constant with height, labeled C, or varies according to 
equation labeled Z. For the constant Am simula- 

tions, the number following the C is the value of Am. For 
the spatially varying Am calculations, the number after- 
wards (along with the AU) describes the radial location 
of the shearing box in our protoplanetary disk model in 
units of AU. An S (L) following the Am value corresponds 
to a domain size of 2H x AH x 8H {8H x 16H x 8iJ), 
which is smaller (larger) than the 4iJ x 8H x 8H size 
of most of the constant Am calculations. The "starter" 
simulation for the AH x 8H x 8H runs is also included in 
the table, labeled Cle5. Finally, all of our calculations 
are carried out at a resolution of 36 grid zones per H . 

3. RESULTS 

3.1. Constant Am Calculations 

We begin by applying some standard diagnostics to 
the set of calculations with constant values of Am. The 
first such diagnostic is the density- wei ghted Maxwell and 
Reynolds stresses (see equation (37) of lBalbus fc HawlevI 
[19981) ■ defined as 



Table 1 

Shearing Box Simulations 



{pVxSVy 



Brr,B,i 



(p) 



(17) 



where the angled brackets denote a volume average over 
the whole domain. Figure [2] shows the time evolution of 
this total stress for the 4H x 8H x 8H runs, normalized 
by the square of the sound speed. The dashed line indi- 
cates the averaged value (from orbit 25 to 53) of the run 
with Am = 10^ from which all of the other 4iJ x8Hx 8H 
runs were restarted. The different Am values are denoted 
by the color. The runs with Am > 1 appear to adjust 
on a roughly 50 orbit timescale after which a statisti- 
cal steady state follows. In general, the stress increases 
with increasing Am (decreasing diffusion) for these runs. 
However, the Am = 10 and Am = 100 runs have roughly 
the same values at late times, as do the Am = 300, 1000, 
and 10000 runs. 

The Am = 1 case deserves extra attention. From 
Fig. [21 it would appear that the turbulence completely 
dies away. A closer examination of the stress histo- 
ries show that the Maxwell stress levels out to a small, 
but positive value, while continuing to slowly decrease 
with time. The Reynolds stress approaches an oscilla- 
tory behavior which occasionally brings it below zero. 
Space-time plots of various quantities in this run indi- 
cate that the gas is no longer MRI turbulent. The rem- 
nant Maxwell stress is the result of a residual large scale 
Bx and By field near the mid-plane, and the Reynolds 
stress appears to arise from residual waves propagating 
through the box. The longer term behavior of this Am 

^ We choose a different restart time for these calculations because 
we decided midway through running our simulations that a larger 
number of cores is significantly more efficient for the variable Am 
runs. Thus, we had to redo the Am = 10^ calculations, and orbit 
22 was chosen because the density-weighted stress at this time 
was roughly equal to orbit 50 of the lower core version of this 
calculation. 
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Ambipolar Diffusion 


Domain Size 


Q 






{Lx X Ly X Lz)H 




CI 


Am = 1, constant 


4x8x8 


decayed 


CIO 


Am = 10, constant 


4x8x8 


0.0046 


ClOL 


Am = 10, constant 


8 X 16 X 8 


0.0055 


ClOOS 


Am = 100, constant 


2x4x8 


0.070 


ClOO 


Am = 100, constant 


4x8x8 


0.0062 


C300 


Am = 300, constant 


4x8x8 


0.024 


ClOOO 


Am = 1000, constant 


4x8x8 


0.022 


ClOOOO 


Am = 10^, constant 


4x8x8 


0.025 


Cle5 


Am = 10^, constant 


4x8x8 


0.038 


Z30AU 


Am at _R = 30AU 


8 X 16 x 8 


0.0016 


ZIOOAU 


Am at ij = lOOAU 


8 X 16 X 8 


0.0015 



value could not be examined because even with STS, 
the diffusion limited time step is very small for Am = 
1; running it out further would be very computation- 
ally expensive. However, these results strongly indicate 
that the MRI turbulence has comp letely decayed away, 
consistent with the results of Bai fc Stond ()2011| ). This 
behavior will play an important role in the variable Am 
simulations of Section 13.21 

We time-average this normalized stress and define the 
Shakura-Sunyaev a parameter. 



(18) 



where the overbar denotes the time average, which is 
done from orbit 100 onwards for most of the constant 
Am runs with Am > 1; from orbit 72 onwards for runs 
ClOL, Z30AU, and ZIOOAU; and from orbit 25 to 53 
for the Am — 10^ "starter" simulation. Fig. [3l displays 
a versus Am for these runs. The arrow on the Am = 
1 run indicates that that the stress level is continually 
decreasing. The trend of a with Am can be compared 
with the u nstratified simulations shown in Fig. 10 of 
iBai fc Ston e (2011) . These tre nds roughly agree, though 
the results from IBai fc Stond ()2011t ) show a monotonic 
increase of a with Am, whereas our results show that 
different Am values can lead to very similar a values. 

This difference may be attributable to different back- 
ground magnetic field strengths as the back ground field 
evolv es via the usual M RI dynamo (e.g., jPavis et al.l 
20M iSimon et al.ll2011b| ). To understand this further, 
first let us consider the space time diagrams of the 
toroidal field By for CIO, ClOO, C300, and ClOOO as 
shown in Fig. [J In these diagrams, the field has been 
averaged over x and y and is plotted in the {t, z) plane. 
The most obvious feature from these diagrams is that 
the period of the dynamo flipping of By changes with 
Am; as ambipolar diffusion becomes stronger, the period 
increases. In particular, for Am = 10, the period is ^ 50 
orbits, and for Am = 100 (only considering times past the 
initial 50 orbit transient period), the period is ~ 15 — 20 
orbits. For Am > 300, the period is ^ 10 orbits as is 
usually observed in stratified MRI simulations. 

The most relevant feature here, however is that the 
background toroidal field strength i s diffe rent for differ- 
ent values of Am. In IBai &: Stond ()2011[ ). it was found 
that with zero net vertical flux, the stress level increases 
with increasing net toroidal flux (which is mostly con- 
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Figure 2. Density- weighted volume average of ttie total (Maxwell 
and Reynolds) stress, normalized by the square of the sound speed, 
versus time for the standard 4H X 8H X SH simulations. The 
magenta line corresponds to Am = 1, blue to Am = 10, red to Am 
= 100, orange to Am = 300, black to Am = 1000, and green to 
Am = 10*. The horizontal dashed line corresponds to the time- 
averaged (from orbit 25 to 53) stress value for Am = 10^ from 
which the other runs were initialized. After an initial transient of 
~ 50 orbits, the simulations with Am> 10 are sustained. There is 
a general trend of increasing stress level with increasing Am. The 
Am = 1 case has turbulence that decays away rapidly. 
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Figure 3. Time-averaged total stress (i.e., a) as a function of 
Am for the standard 4H X 8H X SH simulations. There is a general 
trend of increasing stress level with increasing Am. 

served in unstratified simulations). Therefore, in our 
stratified simulations, two effects are expected to deter- 
mine the saturated stress values: the value of Am and the 
background toroidal field strength, set by the dynamo. 
To demonstrate this effect more robustly, we calculate a 
version of the plasma /3 for the background toroidal field, 



Py EE 2{P)/{Byy^ 



(19) 



where the overbars indicate a time average (from orbit 
100 onwards) and the angled brackets denote an average 
over X and y. This quantity is representative of the am- 
plitude of the oscillating background toroidal field. is 
a function of z only, and we plot it along with the vertical 



profile of the total stress (which has again, been averaged 
in time and for all x and y) in Fig. [5l 

The stress profile reveals the same behavior as that 
in Fig. |3l there is a general trend of increasing stress 
with increasing Am. Furthermore, this increase occurs 
uniformly across all z. However, C300 and ClOOO have 
roughly the same stress profiles, and CIO peaks at around 
the same value as ClOO. Examining the (3y for these par- 
ticular simulations, we see that C300 has a lower value 
(stronger field) than does ClOOO. Similarly, CIO has a 
significantly smaller Py than ClOO. These results con- 
firm that it is indeed the larger background toroidal field 
strength that make the stress levels in run CIO approach 
that in run ClOO, and the stress in run C300 approach 
that in run ClOOO. We note, however, that ClOOO and 
ClOOOO have both the same (3y profiles and the same 
stress profiles. This could indicate that for Am > 1000, 
the turbulence levels are approaching that of ideal MHD. 
The slightly higher stress for Cle5 would then be ex- 
plained by its lower /3j,. 

It remains unclear why the background field strength 
varies in the way that it does. Could this also be con- 
trolled by the value of Am? This is not unreasonable 
considering that ambipolar diffusion already affects the 
period of the toroidal field flipping. The question of ex- 
actly how Am and the dynamo are related is very open. 
Unfortunately, exploring it in detail would take us too 
far from our goals in this paper, and so, we leave it for 
future work. 

The final diagnostic we employ is the two-point au- 
tocorrelation function first used in the context of MRI 
simulations in lGuan et al.l ()2009[ ). We employ th is diag- 
nostic for similar reasons as those in lSimon et al.l ([2012); 
we wish to determine the structure of the turbulent mag- 
netic field and check that the domain sizes we use are 
sufficiently large to properly capture important turbu- 
lent scales. Thus, we define the autocorrelation function 
of the i'^ component of the perturbed magnetic field as 



ACF{SB,{Ax)) 



J SBi{t,x)SBi{t,x + Ax)d^x 
JlB~{t^xf¥x ' 



(20) 



where 5Bi is the value of Bi after subtracting off the 
horizontal mean field. In equation form, 



5Bi{x,y,z) = Bi{x,y,z) - {Bi)xy{z), 



(21) 



and the average denoted by Qxy is the horizontal average. 
Note that we have defined the ACF to be normalized by 
its maximum value (at Ax = Ay = Az — 0). The ACF 
of the total turbulent magnetic field is then defined as 
ACF{SB) = ACF{5Bx) + ACF{SBy) + ACF{SB,). The 
overbar denotes a time average done from orbit 100 to 
125 in aU cases. 

From the figure, it appears that the Am = 100 ACF is 
roughly consistent between a domain size of 2H x AH x 
8H and AH x 8H x 8H, though there is a slight difference 
in the size of the tilted centroid. However, as we will see 
shortly, 2H x AH x 8H is actually too small for Am = 
100. The standard box size, AH x 8H x 8H, appears 
to properly contain the ACF for Am = 100, but not as 
well for Am = 10. The centroid of the Am = 10 case is 
larger and appears to have a longer tail that intersects 
the toroidal boundary. Going to an even larger domain, 
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Figure 4. Space-time plot of By averaged over x and y for Am = 10 (upper left), Am = 100 (upper right), Am = 300 (lower left), and 
Am = 1000 (lower right). The "butterfly" dynamo is present in all simulations, but the period of the By flipping increases with decreasing 
Am. In particular, the period is ~ 40 — 50 orbits for Am = 10 and ~ 15 — 20 orbits for Am = 100. For the other two cases, the period is 
~ 10 orbits, equal to that in ideal MHD calculations. 



8Hxl6Hx8H,the ACF(B) for Am = 10 fooks more well 
contained, though the very end of the tail does appear 
to touch the toroidal boundary. This effect is not as dra- 
matic as in the smaller domain. Going to an even larger 
domain and running Am = 10 is prohibitively expensive 
given our current computational resources. 

Returning to the smallest domain, we note some odd 
features. Despite the reasonable ACF, an inspection of 
the stress history and the a value (see Table [T]) show 
this calculation to be quite different than Am = 100 at 
larger domain sizes. An examination of the space-time 
evolution of By, Fig. [71 brings the point home further, 
as it indicates that the dynamo behavior is completely 
shut off for this particular run. This is again inconsistent 
with the larger domain Am = 100. Thus, we conclude 
that 2H X AH x 8H is too small of a domain for Am = 
100 and will likely be too small for smaller values of Am 
as well. 

As indicated by these results, the standard box size of 
AH X 8H X 8H is large enough for all of our simulations 
except for Am < 10. It is computationally too expen- 
sive for us to run all of our simulations at the larger 
8H X 16H X 8H. So, we elect to use the smaller size 
as our standard, and we compare the evolution of the 
stress between the AH x 8H x 8H and 8H x 16H x 8H 
domains for Am = 10 to justify using a smaller domain 
for comparison between Am = 10 and other Am values. 



Figure [8] shows the Wi?^ evolution for these two domain 
sizes with Am = 10. The use of a smaller domain size 
does not appear to make a difference for this value of 
Am. Furthermore, the By space-time plot for the larger 
domain looks very similar to the smaller domain. Evi- 
dently, we can get away with using a smaller domain for 
Am = 10. However, these ACFs suggest that caution 
be used when running ambipolar diffusion shearing box 
calculations. 

One final thing to note from ACF(i3) is that as Am is 
decreased, the tilted centroid component appears to be- 
come more elongated (hence the need for larger domains) 
and less tilted with respect to the y axis. 

To summarize this section, we find the turbulent stress 
level dependence on Am in vertically stratified simula- 
tions to be generally c onsistent with the results of un- 
stratified simulations (|Bai fc Stond 1201 ID ; a increases 
with Am, and for Am < 1, there is no turbulence. We 
also find that as Am is decreased, larger domain sizes 
are needed in order to properly capture the turbulent 
structures represented by the ACF. These runs serve as 
a baseline for interpreting the results from the next sec- 
tion, in which Am varies spatially and temporally based 
upon our model for surface layer ionization. 



3.2. Variable Am Calculations 
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Figure 5. Vertical profiles of /3y as defined by equation II19I I (top) 
and the total stress normalized by the square of the sound speed 
(bottom). The quantities have been averaged over x and y and 
over time from orbit 100 onwards as described in the text. The 
blue line corresponds to Am = 10, red to Am = 100, orange to Am 
= 300, black to Am = 1000, green to Am = lO**, and dashed line 
to Am = 10^ (the time average for this run is done from orbit 25 to 
the end of the calculation). The stress increases with Am roughly 
uniformly at all heights. There is no obvious trend between fiy and 
Am. 



We now turn our attention to the two calculations with 
spatially and temporally varying Am (the "Z" simula- 
tions in Table [J). As stated previously, these simula- 
tions adopt more realistic non-constant Am values and 
directly address the questions such as "how vigorous is 
outer disk MRI-driven turbulence?" or "what is the mass 
accretion rate in the outer disk?" under the assumption 
that the outer disk is not threaded by any net vertical 
magnetic field. We run them all at the largest domain 
size, 8Hxl6Hx8H, because there are regions of Am < 1 
in these calculations. 

As before, we begin by examining the density- weighted 
stress normalized by the square of the sound speed, as 
shown in Fig. |9l It is clear from the figure that the 
stress levels are lowered dramatically compared to the 
ideal MHD case. Indeed, by calculating the average of 
the curve from orbit 72 onward, we find that the a val- 
ues are ~ 10""^, one order of magn itude below wh at is 
expected from observations (Hartm ann et al.lll998l ). 

While the turbulence levels initially decrease drasti- 
cally, it has not completely died away. Consider the 
space-time diagrams in the {t,z) plane for various hor- 
izontally averaged quantities of run ZIOOAU, as shown 
in Fig. [TOl From both the Maxwell and Reynolds stress 



plots, there is a significant region in z over which the MRI 
is indeed active. The vertical structure is consistent with 
what we would expect from our ionization profile; there 
is a significant "ambipolar dead zone" corresponding to 
where Am = 1 and the higher Am regions correspond 
to turbulent activity. We will calculate an actual mass 
accretion rate at both 30 AU and 100 AU in Section HTTl 

Despite there being very little Maxwell stress in the 
ambipolar dead zone, there is still Reynolds stress in this 
region. This stress is likely produced by the active layers, 
similar to what is o bserved in simulations th at include an 
Ohmic dead zone ([Fleming fc Stonel[200l . It IS possi- 
ble, however, that as in run CI, the Reynolds stress here 
is simply left over from the turbulent state from which 
the run was restarted. To test this notion, we restarted 
ZIOOAU at orbit 85 and set all perturbed velocity com- 
ponents to zero throughout the domain (the shear profile 
is of course maintained). We find that the active MRI 
layers do indeed induce velocity fluctuations within the 
ambipolar dead zone, which leads to a positive Reynolds 
stress on average. 

Finally, the horizontally averaged By space-time plot 
shows interesting behavior. Within the ambipolar dead 
zone, the toroidal field remains fixed in sign, though the 
magnitude appears to be decreasing. Between 1 and 2 
scale heights above and below the mid-plane, the usual 
MRI dynamo reappears, with a By oscillation period of 
^ 10 orbits, identical to ideal MHD. Again, from the 
Maxwell stress plot, it is obvious that this region is tur- 
bulent. For \z\ > 2H, however, there is a strong toroidal 
field that remains stationary. There are slight changes in 
magnitude as the toroidal field from the turbulent region 
propagates outwards. 

Examining the same diagrams for Z30AU reveals very 
similar behavior. There is a thin layer of strong Maxwell 
and Reynolds stress located around \z\ ^ 2H, along with 
a positive Reynolds stress induced in the ambipolar dead 
zone by the active layers. In this calculation, however, 
the active region is much narrower in z when compared 
to ZIOOAU. This is because the FUV photons do not 
ionize the disk quite as deep as at i? = 100 AU, and as 
is usual, the MRI appears to be inactive for \z\ > 2H . 
Thus, the MRI- active region is forced to be contained 
within a smaller z region. 

One thing t o question is whethe r or not the phe- 
nomenology of I Simon et al.l (|2011bl ) could play a role 
here if we were to integrate these simulations further. 
In that paper, the authors found that the toroidal and 
radial field in a resistivity dominated mid-plane region 
would evolve in time. Occasionally, the toroidal field 
would grow to strong enough levels that the MRI would 
be temporarily reactivated, after which the Ohmic re- 
sistivity would quench the turbulence again. In these 
simulations, we do see that the toroidal field in the am- 
bipolar dead zone changes in amplitude over time. Is it 
possible, then, that the field could grow to large enough 
values to eventually re-activate the MRI in that region? 
Integrating our simulations out for many hundreds of or- 
bits is prohibitively expensive, and therefore, we will not 
be able to observe any such variability in our stratified 
simulations. Instead, we have run an unstratified shear- 
ing box of size 8H x 16H x H with Am = 1 and uniform 
toroidal field of strength /3 = 10, chosen to determine if 
a strong toroidal field MRI will be active with Am = 1. 
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Figure 6. Autocorrelation function (ACF) of tlie magnetic field, as defined by equation II20I I. for simulations (from left to right) with Am 
= 100 and size 2HxiHx 8H, Am = 100 and size 4H x8H X 8H, Am = 10 and size 4H x8H X 8H, and Am = 10 and size 8H X 16H X 8H. 
As Am is decreased, larger and larger domains are needed to properly contain the ACF. Furthermore, the tilted centroid feature becomes 
less tilted with respect to the y axis and more elongated as Am is decreased. 
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Figure 7. Space-time plot of By averaged over x and y for the 
Am = 100 run at domain size 2H X AH X 8H. There is a remnant 
dynamo behavior that rapidly dies away as the simulation adjusts 
to the new value of Am. Eventually, the dynamo activity ceases 
altogether, which is inconsistent with the larger domain counter- 
part of Am = 100. This domain size appears to be too small to 
properly capture the dynamo at Am = 100. 

We start the simulation with a relatively large amplitude 
perturbation to the density and velocities, such that the 
initial conditions should already be reasonably nonlinear. 
As a control, we also ran this identical setup with Am = 
10^. With Am = 1, we observed decay of the initial per- 
turbations and no development of any MRI turbulence, 
whereas with Am — 10^, the MRI becomes active and 
generates sustained turbulence for many orbits. Thus, 
even in the presence of a strong field, Am = 1 is suffi- 
cient to quench MRI-driven turbule nce. We therefore do 
not expect the behavior observed in lSimon et al.l ()2011bO 
to occur in these simulations. 

The vertical structure results can also be summarized 
in the time- and horizontally-averaged quantities as a 
function of z, as is shown in Fig. [TT] From the top two 
panels, it is obvious that there is a double-peak struc- 
ture to the stress; no doubt a result of the ambipolar 
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Figure 8. Density-weighted volume average of the total (Maxwell 
and Reynolds) stress, normalized by the square of the sound speed, 
versus time for the Am = 10 simulation at a domain size of 
AH X 8H X 8H (solid line) and 8H X 16H X 8H (dotted line). 
Note that the 8H X 16H X 8H run was restarted from orbit 22 
of a different "starter" simulation. For comparison purposes, we 
translated the solution to the right by 28 orbits. The curves show 
a nearly identical evolution. 

dead zone. Within this dead zone, the Reynolds stress 
dominates over the Maxwell stress. From the space-time 
diagrams above, any nonzero Maxwell stress within this 
region likely results from large scale correlations in the 
magnetic field rather than any sort of turbulence. Note 
that there are regions where the stress can go negative, 
and since the vertical axis is logarithmic, the curves are 
simply cut off where the values drop below zero. 

The bottom two panels show the various energies (i.e., 
thermal, kinetic, and magnetic) as a function of height. 
The thermal pressure dominates over the vast majority 
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Figure 9. Density-weighted volume average of the total (Maxwell 
and Reynolds) stress, normalized by the square of the sound speed, 
versus time for the Z30AU simulation (solid line) and the ZIOOAU 
simulation (dashed line). The curves show a clear drop from the 
initial state of vigorous turbulence to levels near Wn^/c'^ 10^''. 

of the disk's vertical structure. However, for \z\ > 2.5H, 
the magnetic energy dominates over all other energies. 
We note in particular the very strong magnetic domi- 
nance in the upper regions of Z30AU (lower left panel). 
This magnetic field is stationary in time according to the 
space-time diagram for this run, akin to the lower panel 
of Fig. 1101 a net toroidal field remains stationary in time 
for \z\ > 2H. 

We attribute some of this behavior to the gas den- 
sity floor. Indeed, looking at the vertical pressure profile 
(which follows the gas density) , the pressure is prevented 
from going below 10^^ of its initial peak value. The re- 
sulting gradient in the gas pressure has an effect on the 
buoyant properties of the field. We address this issue 
further in the next section. 

3.3. The Effect of the Density Floor 

To test the effect of lowering the density fioor, we 
restarted simulations CIO and Z30AU at orbits 120 and 
56 respectively and lowered the density floor to 3 x 10~^, 
roughly 3 times smaller than the original density floor. 
While this new value for the floor may still be higher 
than what the density would naturally be, it becomes ex- 
tremely expensive to run lower density floor calculations 
even for a short integration. Thus, our main goal in run- 
ning these calculations is to observe the immediate effect 
of lowering the floor on properties such as the density- 
weighted stress and the buoyancy of the magnetic field 
in the upper disk regions. Is the evolution of the system 
altered drastically? What, if any, changes occur? 

Examining the space-time diagrams for these restarted 
runs indicate that the lowered density floors lead to en- 
hanced buoyancy of the magnetic flelds. For example, 
consider the lower left panel of Fig. [11] Once the floor 
is lowered, the magnetic energy for z < —2H immedi- 
ately drops to lower values as field rises away from the 
mid-planeH The same behavior is observed in run CIO 
for both z < —2H and z > 2H . Thus, the presence of 

The field for z > 2H does not appear to change significantly 
over the time that we integrated this simulation. However, we did 



a strong magnetic field for \z\ > 2H, as is shown in the 
lower panel of Fig.[TOl is artificial. That bein g said, based 
upon t he results of ideal MHD calculations ( Simon et al.l 
l2011b( ) in which the density is (on average) larger than 
the fioor value at all locations, we still expect the field 
to be superthermal in this region; it will just be weaker. 

Does this feature affect our main results so far? An ex- 
amination of the stress evolution shows that the decrease 
in the density floor does lead to a decrease in the volume- 
averaged stresses. Quantitatively speaking, in Z30AU, 
the magnetic energy averaged over z < —2H drops by a 
factor of ~ 4 in going from orbit 56 to 78.5. During this 
same time, the volume-averaged stresses decrease by a 
factor of ^ 1.5, and the stress does not appear to be lev- 
eling off. Running the simulation out further is compu- 
tationally prohibitive given the small time step incurred 
by the lower density floor. So, there is an effect due to 
the density floor. Since we have not fully quantified this 
effect, it should be borne in mind when we calculate mass 
accretion rates in Section 14.11 The question of whether 
or not the density floor influences the turbulent velocities 
at large \z\ is addressed in Section 

4. IMPLICATIONS FOR PROTOPLANETARY DISK 
STRUCTURE AND EVOLUTION 

4.1. Mass Accretion Rate 

One of the most important properties of disk evo- 
lution is the mass accretion rate due to angular mo- 
mentum transport. We can calculate this quantity, M, 
for Z30AU and ZIOOAU by utilizing equation (40) in 
iBalbus fc Ha^ievl (|l998f ). which assumes accretion is in 
steady-state. We take this equation in the limit that 
R ^ Ro, slS appropriate for the shearing box approxima- 
tion: 

M = '^W^ ~ 8.5 X lO-^ai?^^' Mg/yr, (22) 

where the expression on the right comes from applying 
this formula to the minimum mass solar nebula model, 
and a i s deflned as in equation (jl8p . The deflnition of 
Wr^ in IBalbus fc HawlevI (|1998D . is the same as ours m 
the sense that it is a density weighted stress. From this 
equation, we calculate that M « 2.5 x 10^^ Mq /yi at 30 
AU and M w 1.2x IQ-^Ma/yr at 100 AU. Of course, due 
to the effect from the high density floor (see discussion 
in Section [231) , the actual accretion rates are likely to be 
lower, and these values serve as an upper limit. 

These accretion rates are at least one order of magni- 
tude too small (likely even smaller, again due to the den- 
sity floor) to account for the ob served accretion rates in 
Classical T-Tauri systems (e.g.. iHartmann et al.lllQQSh . 

We can also con ipar e our res ults to semi-analytical pre- 
dictions made bv iBail ()2011al) based upon the results of 
unstratiflcd simulations with a net vertical magnetic flux 
fBai & Stone 2011). This comparison will allow us to 
gauge the potential importance of having a net vertical 
flux in regions of strong ambipolar diffusion. By flrst ex- 
tracting the vertical proflles of p and Am, and then by as- 
suming constant magnetic field strength across the MRI 
active region of the disk, we can estimate the accretion 

not evolve this simulation for very long, and we speculate that a 
longer evolution would show a change in the magnetic field strength 
in this region. 
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Figure 10. Space-time plots of Maxwell stress (top), Reynolds stress (middle) and toroidal field, By (bottom) for run ZIOOAU, where 
each quantity has been averaged over all x and y. The stresses are normalized by the initial peak gas pressure (Po = 5 X 10"''), whereas By 
is in code units. There is a clear "active" region between 1 and 2 scale heights above and below the mid-plane, where the Maxwell stress is 
non-negligible and the toroidal field dynamo is active. Within 1 scale height of the mid-plane, there is no turbulent activity, though there 
is non-zero Reynolds stress. 
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Figure 11. Time- and iiorizontally-averaged stresses (upper row) and energies (bottom row) vs. z in scale heights for the run at 30 AU 
(left column) and 100 AU (right column). The time average is done from orbit 72 onwards, and the horizontal average is done over all x 
and y. In the stress plots, the solid line is the Maxwell stress, and the dashed line is the Reynolds stress. In the energy plots, the solid line 
is the gas pressure, the dashed line is the magnetic energy, and the dotted line is the kinetic energy. All quantities are normalized by the 
initial peak gas pressure (Po = 5 x 10"'^). 
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rate for any given field strength using equation (28) of 
[Bail (|2011aD . From this approach, M = 9.8 x lO^^Mg/yr 
and 3.5 x lO-^Mg/yr at 30 AU and 100 AU, respectively, 
roughly one order of magnitude larger than our calcu- 
lated rates and in general agreement with observations. 

These estimates suggest that it is very likely that a 
net vertical field is required to attain the necessary tur- 
bulence levels, if the MRI is indeed the dominant mech- 
anism by which angular momentum is transported. We 
will carry out actual shearing box simulations with verti- 
cal stratification and ambipolar diffusion in the presence 
of a net vertical magnetic field in Paper II. 

4.2. Turbulent Linewidth 

Another property of disks of recent interest has been 
the density-weighted distribution of turbulent velocities 
as a function of disk radius and h eight above the mid- 
plane. This was first calculated bv iSimon et al.l ()2011aD 
for local MRI simulations without ambipolar diffusion, 
but in cluding the e f fects of Ohmic resistivity. Another 
study, iForgan et al.l (|2012l ) , carried out a similar analysis 
for global calculations of self-gravity driven turbulence. 
These distributions are a first order approach to making 
a connection with observational constraints of turbulent 
line broadening, suc h as was observed in the sub-mm by 
iHughes et al.l pOllD . In particular, the density- weighted 
velocity represents the probability of observing a line 
with a particular turbulent velocity broadening. 

Here, we car r y out a n identical analysis to that done 
in ISimon et al.l ()2011af) for both of our variable Am cal- 
culations (see that paper for the exact details of how 
to calculate the velocity distribution). Figure fT2l shows 
this velocity distribution for Z30AU (top and middle) 
and ZI OOAU (bottom). As was observed in ISimon et al.l 
(|2011al ) for their calculations with a strong Ohmic dead 
zone (see their Fig. 4, top panel), we also observe a rea- 
sonably strong gradient in peak velocity as one probes 
deeper towards the mid-plane. Indeed, the mid-plane 
velocity distribution peaks around v/cg ~ 0.01, just as in 
the Ohmic case. Furthermore, as one probes the surface 
layers of the disk \z\ > 3H, the peak of the distribution 
occurs around v/cs ^ 0.2 — 0.4. There is also a non- 
negligible supersonic tail; in ZIOOAU, this component 
comprises ~ 1% of the horizontal and vertical distribu- 
tions at z > 3H. In Z30AU, this component comprises 
^ 2% of the vertical distribution and ^ 7% of the hori- 
zontal distribution at z < — 3-ff. 

One thing to note from the top panel of the figure is 
that the red and black curves nearly lie on top of each 
other. This is likely an artificial effect resulting from 
the very dominant magnetic field that is stationary for 
z > 2H in Z30AU (which itself results from the relatively 
large density fioor applied in this calculation) . The mag- 
netic field is not nearly as dominant for z < ~2H (the 
other side of the mid-plane), and so the middle panel of 
Fig. [12] shows the velocity distribution calculated from 
this side. The distribution looks much more similar to 
the other distributions. 

To further test the effect of the density fioor on our 
velocity distributions, we have rerun Z30AU with the 
density fioor lowered to 3 x (as discussed in Sec- 

tioning!]). We calculated the velocity distribution for this 
region during two separate periods, each averaged over 8 



orbits. We do not see any significant difference between 
these velocity distributions and that shown in the middle 
panel of Fig. [12] 

Finally, the rough agreement between turbulent veloc- 
ity distributions for the vertical velocity and the "in- 
plane" velocities suggest that turbulent mo tions will be 
isotrop ic, consistent with previous results ([Simon et al.l 
l2011al) . We point out that velocity distribution for z > 
in ZIOOAU suggests that the fiow is slightly anisotropic. 
We are not entirely sure why this is the case. However, 
when we restarted this run and set the turbulent veloc- 
ity to zero (as described in Section 13. 2|) . the resulting 
velocity distribution was again nearly isotropic. Since 
this isotropy is present in all of the other cases, it seems 
more likely that the distribution for z > in ZIOOAU is 
a peculiar case, perhaps resulting from the exact nature 
of the turbulent state from which this run was initiated. 

5. SUMMARY AND CONCLUSIONS 

We have run local shearing box simulations of MRI- 
driven turbulence in the presence of ambipolar diffusion 
and in the absence of a net vertical magnetic fiux. These 
simulations were designed to address two primary ques- 
tions: 

• How does MRI-driven turbulence behave in the 
presence of both ambipolar diffusion and vertical 
gravity? 

• What are the implications for turbulence in the 
outer regions of protoplanetary disks where am- 
bipolar diffusion is dominant? 

With the ambipolar Elsasser number. Am, (see equa- 
tion [5]) remaining constant, we addressed the first ques- 
tion. We found that the density-weighted stress de- 
creases with increasing ambipolar diffusion and decays 
to negligible values for Am < 1. Another parameter that 
controls the stress levels, however, is the amplitude of 
the toroidal field strengt h, which oscillates in time due to 
the dynamo (jSimon et a l. 2011b). Namely, the stronger 
the amplitude of oscillation (but not too strong which 
would suppress the MRI in the presence of ambipolar 
diffusion), the larger the turbulent stresses. The average 
strength of this varying field combined with ambipolar 
diffusion affect the MRI in such a way that the turbu- 
lent stress does not necessarily increase monotonically 
with decreasing diffusion. All our results a re consistent 
with unstratified numerical simulations of iBai fc Stond 
(|2011f ). although the subtleties with the ambipolar MRI 
dynamo due to the addition of vertical gravity does not 
guarantee a one-to-one correspondence between the level 
of diffusion and that of the turbulent stress. 

Additional noticeable effects emerged from these con- 
stant Am calculations. First, as ambipolar diffusion is 
increased, the dynamo oscillation period becomes longer. 
Above Am ~ 100, the oscillation period approaches the 
ideal MHD limit of 10 orbits. This result opens up more 
questions than it answers, as we do not yet have an un- 
derstanding of the dynamo in the ideal MHD limit, let 
alone including a diffusion term. However, it may be 
insightful to apply an ambipolar diffusion term to cur- 
rent, simplified models of the MRI dynamo to help in 
further understanding the physics of the dynamo. Sec- 
ond, the typical turbulent fiuctuations become larger in 



15 



R = 30 AU 




0.001 0.010 0.100 1.000 
R = 30 AU 



10.000 




0.001 0.010 0.100 1.000 
R = 100 AU 



10.000 




0.001 0.010 0.100 



1.000 



10.000 



Figure 12. Density-weighted turbulent velocity distributions for 
Z30AU (top and middle) and ZIOOAU (bottom). The top and 
bottom panels correspond only to velocities at z > 0, whereas 
the middle panel corresponds to z < 0. Each color corresponds 
to different depths over which the distribution is calculated, as 
labeled. The dashed lines are the vertical turbulent velocity \vz\/cb, 
and the solid lines are the azimuthally averaged disk planar velocity 
Dftl/cs. Relatively small velocities exist towards the mid-plane, as 
turbulence is very weak there. However, far from the mid-plane, 
the velocity consistently peaks at v/cs ~ 0.2 — 0.4, with part of the 
distribution going into the supersonic regime. 



scale (i.e., have a larger correlation length) and become 
more aligned with the azimuthal direction as ambipo- 
lar diffusion is increased. This last point has important 
implications for local simulations that include ambipo- 
lar diffusion. Increased diffusion evidently favors larger 
scale fluctuations, and for simulations with strong diffu- 
sion, larger shearing boxes are required. A box size of at 
least 8H x 16H x 8H is required to properly capture, the 
MRI turbulence with Am < 10. This also motivates fur- 
ther studies of ambipolar diffusion and the MRI in global 
simulations, where one is not limited by scales of order 
H. 

To answer the second of our motivating questions, we 
ran additional simulations that included a physically mo- 
tivated model for the ionization structure (and hence Am 
profile) of the protoplanetary disk. These runs include 
the effect of a strong FUV ionizati on laye r based upon the 
work of Perez-Becker &: Chiang ()2011bl ). where a large 
ionization fraction / ~ 10~^ can be achieved in a very 
thin layer above and below the disk mid-plane, while the 
rest of the disk is assumed to have Am = 1 in accordance 
with the chemistry calculations of jBai! (12 01 la [) . Although 
this ionization model still bares uncertainties, it provides 
the essential physical ingredients that allow us to explore 
the gas dynamics in the outer regions of protoplanetary 
disks in a realistic manner. 

We find that despite this strong FUV ionization, the 
mass accretion rate is of order lO~^M0/yr, too small to 
accou nt for observed accretion rates measured in T T auri 
stars (|Gullbring et all [19981 iHartmann et al.lll998f) . In 
fact, this estimate should be treated as an upper limit due 
to the increase in stress from the relatively large density 
floor employed. This small accretion rate and the pres- 
ence of the ambipolar dea d zone is remin iscent of models 
for the Ohmic dead zone (|Gammi i [l99l . many of which 
also yield low accretion rates. The problem posed by an 
ambipolar dead zone is, however, more serious. Close 
to the star, the viscous time scale R^/v can be short 
compared to the disk lifetime given even a weak residual 
stress. If an Ohmic dead zone can be supplied with gas 
from further out, it is then possible to imagine gas ac- 
cumulating there until some additional i nstability allows 
gas to flow through on to the star (e.g.. lArmitage et al.l 
I200H iZhu et all l2010t iMartin et aLll2012| ). In the am- 
bipolar dead zone, conversely, low levels of stress imply 
that the main mass reservoir is permanently inactive. 
Since this appears to contradict observations, our simu- 
lations are either missing some important aspect of the 
physics, or our basic understanding of disks is incorrect. 

There are several possibilities. First, we could be miss- 
ing some effects that strongly influence the MRI. In- 
cluding the Hall effect may substantially enhance the 
strength of MRI turbulence, as indicated by previous 
w orks (see the In t roduct ion). Except for the simulat ions 
of ISano fc Sto^ (|2002aD and lSano fc Stmi (|2002bD . no 
one has carried out non-linear MRI simulations of the 
Hall effect. Somewhat discouragingly, though, the re- 
sults of the Sano studies suggest that the Hall effect 
does not have a strong influence when Ohmic resistiv- 
ity is dominant. The same could conceivably be true for 
ambipolar diffusion dominated regions of the disk as this 
diffusion acts similar to Ohmic resistivity; it damps out 
the turbulence. However, this most definitely needs to 
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be investigated in future work. 

Alternatively, it is possible that angular momentum 
transport in the outer disk is dominated by an entirely 
different physical mechanism, though the known candi- 
dates are not expected to be efficient in this region unless 
the di sk mass is high enough that self-gravity is signif- 
icant ([Armitagd 1201 1[ ). It could also be the case that 
the surface density drops off less rapidly than that as- 
sumed in our model. If this is the case, then the induced 
Reynolds stress from the active layers may be able to 
transport more angular momentum outwards, thus in- 
creasing the accretion rate. However, a larger column in 
the outer disk leads to a smaller active region and thus a 
smaller induced Reynolds stress. Finally, one cannot ex- 
clude the possibility that the disk on 30-100 AU scales is 
genuinely inviscid, with observed accretion coming from 
a larger-than-expected reservoir of gas closer to the star. 

However, the most promising possibility for explaining 
these low accretion rates is the exclusion of a net vertical 
magnetic flux, which is the most optimal m agnetic geom- 
etry for the MRI with ambipolar diffusion ()Bai fc Stond 
1201 If ). Net vertical field allows the MRI to operate at ar- 
bitrarily small values of Am and permits relatively strong 
turbulence with a ^ .01 at the mid-plane region of the 
outer disk ()Baill201l"airbl ). It should also make the accre- 
tion in the FUV ionized layer much more efficient since 
the vertical field will be relatively strong here due to the 
drop in gas pressure away from the mid-plane. Indeed, 
an estimate of the a ccretion rate ba sed on the models o f 
iPerez-Becker fc Chiang. (2011b) andlBai fc Stond (poTTI ) 
and the simulations o f iBai fc Stond (I2011D (which include 
the effects of a net vertical magnetic flux) returns a much 
more optimistic M ^ lO~^M0/yr. These considerations 
strongly motivate our companion paper, where we will 
study in detail the effect of net vertical fields on the ac- 
cretion process in the limit of strong ambipolar diffusion. 

The outer regions of protoplanetary disks can be re- 
solved at sub-mm wavelengths, and with this in mind we 
have calculated the probability distribution for turbulent 
velocities as a function of height within the disk. We find 
that, although the ambipolar dead zone severely restricts 
the accretion rate, turbulence in the active surface lay- 
ers remains strong. We obtain peak values of ~ 0.4 of 
the isothermal sound speed. These results are similar 
to those found by Simon et al. (2011a) in calculations of 
Ohmic dead zones at smaller radii, and are consistent 
with o bservations of HD 162396 made by iHughes et al.l 
(|2011h . We predict that the turbulent velocity (and 
hence line width) ought to be a strong function of height 
at radii where an ambipolar dead zone is present, which 
may be testable given observations of multiple molecular 
tracers that probe different depths within the disk. 

The most promising avenue for observational progress 
lies in ALMA measurements of the spatial structure 
and velocity field of protoplanetary disks on the same 
(large) scales as those considered theoretically in this pa- 
per. As we have noted, imp roved measurements of tur- 
bulent line broadening (|Hughes et al.ii2()Tll ) at different 
depths within the disk can potentially provide direct con- 
straints o theoretical models. Such observations appear 
to be technically feasible (Hughes, private communica- 
tion). Our results, however, also motivate consideration 
of disk evolution scenarios that are substantially different 



from those usually adopted in the interpretation of ob- 
servational data. It is commonly assumed, for example, 
that the outer edges of protoplanetary disks expand sig- 
nificantly as the disk evolves viscously. If this is true, 
then measurements of the disk surface density profile 
(when fit, for example, by similarity solutions) constrain 
the radial variation of the angular momentum transport. 
Our results, on the other hand, suggest that some disks 
(those with negligible net vertical fields) may not evolve 
viscously at all on large scales. It may therefore be use- 
ful to consider models in which qualitatively different 
pathways of disk evolution are driven by variations in 
the initial magnetic flux, rather than by differences in 
the initial mass and angular momentum content of disk 
gas. Since strong vertical fields can also le ad to angu- 
lar m omentum loss via disk winds (e.g., Sai meron et al.l 
1201 Ih . one possible scenario is that strong vertical fields 
lead to wind-dominated disks, weaker fields to stimulated 
MRI-driven evolution, and no vertical field to effectively 
inviscid disks. 
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